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Abstract 

Transient spike trains consisting of M (= 1 - 5) pulses generated by single Hodgkin-Huxley 
(HH) neurons, have been analyzed by using both the continuous and discrete wavelet transfor- 
mations (WT). We have studied effects of variations in the interspike intervals (ISI) of input 
spikes and effects of random noises on the energy distribution and the wavelet entropy, which 
are expressed in terms of the WT expansion coefficients. The results obtained by the WT are 
discussed in connection with those obtained by the Fourier transformation. 
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1 Introduction 

During the last half century, extensive experimental and theoretical studies have been 
made on functions of brain where neurons communicate information by action potentials 
or spikes. Because of the complexity of spike signals, little is known how information is 
carried by spikes at the moment 0-fl. It has been widely believed that information is 
encoded in the average firing rate of individual neurons (rate code). Andrian M first noted 
the relationship between neural firing rate and stimulus intensity, which forms the basis 
of the rate code. Actually firing activities of motor and sensory neurons are reported to 
vary in response to applied stimuli. In recent years, however, it has been proposed that 
detailed spike timings play an important role in information transmission (temporal code): 
information is assumed to be encoded in interspike intervals (ISIs) or in relative timings 
between firing times of spikes [^]-[|TU|. Indeed, experimental evidences have accumulated 



in the last several years, indicating a use of the temporal coding in neural systems [11 



P~5|1 . Human visual systems, for example, have shown to classify patterns within 250 ms 



despite the fact that at least ten synaptic stages are involved from retina to the temporal 
brain fli~5 |. The transmission times between two successive stages of synaptic transmission 



are suggested to be no more than 10 ms on the average. This period is too short to allow 
rates to be determined accurately. 



1 e- mail : hasegawa® u- gakugei . ac . j p 
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Although much of debates on the nature of the neural code has focused on rate versus 
temporal codes, there are other important issue to consider: information is encoded in 
the activity of single (or very few) neurons or that of a large number of neurons. A major 
question on the population code (or ensemble code) concerns how the information is carried 
by correlation between different action potentials fired by different neurons. The simplest 
assumption is that little or no amount of information is carried by such correlation and 
that information is coded in the relative firing rates of the neurons (ensemble rate code). 
This code model is supported by many experiments and has been adopted in the 
majority of analysis of neural coding |6fl ||18|| . On the contrary, it is assumed that relative 
timings between spikes in ensemble neurons may be used as an encoding mechanism for 



perceptional processing (ensemble temporal code) ||T9|-||2T|. A number of experimental 



data supporting this code have been reported in recent years [22|- ||25|| . For example, data 
have demonstrated that temporally coordinated spikes can systematically signal sensory 
object feature, even in the absence of changes in firing rate of the spikes ||23|| . It is currently 
controversial what kind of code is employed in the real neural systems [|TJ-@- 

This shows the importance of studying how spike signals may convey information with 
the sub-millisecond resolution. The response of spikes to stimulus has been best studied 
when it is periodic in time, for which the Fourier transformation (FT) method is usually 
adopted for its analysis. The FT decomposes a signal into its constituent components. 
The information theory provides a powerful tool for analyzing the nature and quality of 
a neuronal code. A natural approach to quantify the degree of order of a complex signals 
is to consider its spectral entropy, as defined from the FT power spectrum [p6 |. The 



spectrum entropy shows how the FT power spectrum is concentrated or widespread; an 
ordered activity with a narrow peak in the frequency domain yields a low entropy while 
a disordered activity with a widespread frequency distribution leads to higher entropy. 

The FT requires that a signal to be examined is stationary, not giving the time evo- 
lution of the frequency pattern. Actual biological signals are, however, not necessarily 
stationary. It has been reported that neurons in different regions have different firing 
activities. Furthermore even within a given region, firing property depends on its states. 
The typical example is found in thalamus, which is the major gateway for the flow of 
information toward the cerebral cortex. In arousal spikes with gamma oscillations (30-70 
Hz), mainly 40 Hz, are reported, whereas spindle oscillations (7-14 Hz) and slow oscilla- 
tions (1-7 Hz) are found in early sleeping and deepen sleeping states, respectively P7|| . In 
hippocampus, gamma oscillation occurs in vivo, following sharp waves ||28|| . In neo-cortex, 
gamma oscillation is observed under conditions of sensory signal as well as during sleep 
29| . Spike signals in cortical neurons are generally not stationary, rather they are tran- 



sient signals or bursts [ 3(| , whereas periodic spikes are found in systems such as auditory 
sytems of owl and the electrosensory system of electric fish [[32 1. 

The limitation of the FT analysis can be partly resolved by using a short-time Fourier 
transform (STFT). Assuming the signal is quasi- stationary in the narrow time period, 
the FT is applied with time-evolving narrow windows. Then STFT yields the time evo- 
lution of the frequency spectrum. The STFT, however, has a critical limitation violating 
the uncertainty principle, which asserts that if the window is too narrow, the frequency 
resolution will be poor whereas if the window is too wide, the time resolution will be less 
precise. This limitation becomes serious for signals with much transient components, like 
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spike signals. 

The disadvantage of the STFT is overcome in the wavelet transform (WT) |33|]. In 



contrast to the FT, the WT offers the two-dimensional expansion for a time-dependent 
signal with the scale and translation parameters which are interpreted physically as the 
inverse of frequency and time, respectively. As a basis of the WT, we employ the mother 
wavelet which is localized in both frequency and time domain. The WT expansion is 
carried out in terms of a family of wavelets which is made by dilation and translation 
of the mother wavelet. The time evolution of frequency pattern can be followed with an 
optimal time-frequency resolution. 

The WT appears to be an ideal tool for analyzing signals of a non-stationary nature. 



In recent years the WT has been applied to an analysis of biological signals [M\ , such as 
electoencephalographic (EEG) waves |471HM) and spikes [49fl- |p^| . EEG is a reflection 



of the activity of ensembles of neurons producing oscillations. By using the WT, we can 
obtain the time-dependent decomposition of EEG signals to 5 (0.3-3.5 Hz), 9 (3.5-7.5 
Hz), a (7.5-12.5 Hz), j3 (12.5-30.0 Hz) and 7 (30-70 Hz) components g7J-g. It has been 
shown that the WT is a powerful tool to the spike sorting in which coherent signals of a 
single target neuron are extracted from mixture of response signals |HD[-pT|. The WT has 



been probed [||| to be superior than the conventional analysis methods like the principal 
component analysis (PCA) f5"3fl . 

Besides the FT, another useful way to describe the dynamical behavior of stationary 
signals is to use Lyapunov exponents and correlation dimensions |54|. Lyapunov expo- 



nents, which were initially introduced for an analysis of chaos, measure the rate at which 
nearby points on an attractor diverge or converge along nearby trajectories. Correlation 
dimensions are meaningful quantities to examine whether a given signal is chaotic or not. 
We must note, however, that these quantities are defined only for stationary signals. 

It is the purpose of the present paper to make a WT analysis of transient spike signals 
which have not been quantitatively discussed so far. We have analyzed spike trains con- 
sisting of M (=1-5) pulses, generated by the Ho dgkin- Huxley (HH) neuron model ]55| . 
We have calculated the energy distribution and the wavelet entropy which are expressed 
in terms of the WT expansion coefficients. The WT is classified to the continuous wavelet 
transform (CWT) and the discrete wavelet transform (DWT). The former treats the scale 
and translation parameters as continuous variables while the latter adopts the variables 
only at the discrete points. Both the CWT and DWT have advantages. The CWT pro- 
vides us with the intuitively clear results. On the contrary, the DWT, which is based on 
the ortho-normal basis, can be quickly performed by using the multi-resolution analysis. 

Our paper is organized as follows: In the next § 2.1, we describe the HH neuron model 



55fl , which is considered to be the most realistic neuron model among proposed theoretical 



models. In § 2.2. we briefly mention the CWT and DWT, presenting expressions for the 
energy distribution and the wavelet entropy. The results of WT analysis of spike signals 
by using both CWT and DWT are presented in § 3, where the effect of ISI of spike signals 
and of noises are studied. The final S 4 is devoted to discussions and conclusion. 



2 Calculation Method 
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2.1 HH Neuron Model 

In order to generate spike trains, we adopt a single HH model pH|. The neuron is as- 
sumed to be in the environment with the independent Ornstein-Uhlenbeck (OU) noises. 
Dynamics of the membrane potential V of the HH neuron is described by the non-linear 
differential equation given by 



C dV(t)/dt = -I ion (V, m, h, n) + J ps + I n , (1) 



where C — 1 /zF/cm 2 is the capacity of the membrane. The first term of Eq.(l) expresses 
the ion current given by 

I ioa (V, m, h, n) = g^m 3 h(V - K Na ) + g K n\V - V K ) + g L (V - V L ). (2) 

Here the maximum values of conductivities of Na and K channels and leakage are $Na — 
120 mS/cm 2 , g K = 36 mS/cm 2 and g^ = 0.3 mS/cm 2 , respectively; the respective reversal 
potentials are Vn & — 50 mV, Vk = —77 mV and Vl = —54.5 mV. Dynamics of the gating 
variables of Na and K channels, m, h and n, are described by the ordinary differential 
equations, whose details have been given elsewhere []5B |. 

The second term in Eq.(l) denotes the postsynaptic currents given by 

^ ps = 9s {V a - V s ) a(t - t im ), (3) 

with the alpha function a(t): 

a(t) = (t/r s ) e-^ 0(t), (4) 

where 0(t) = 1 for x > and for x < 0. Equation (3) expresses the postsynaptic 
current of the neuron, which is induced by the presynaptic spike-train input given by 

M 

Ui{t) = V & S(t-tim). (5) 

m=l 

where V a means its magnitude, M the number of input pulses and t- im is the m-th firing 
time of the input. The interspike interval (ISI) of input signal is defined by 

Ti m ^im+1 t[ m . (6) 

The third term in Eq.(l) denotes the OU noise given by 

T n dI n /dt = -I n + £(t), (7) 

with the Gaussian white noise 



< m >= o, 



<mm>=2D5(t-t'), (9) 

where the bracket < X > denote the time average. The intensity and the correlation time 
of white noises are given by D and r n , respectively 
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Differential equations given by Eqs.(l)-(9) are solved by the forth-order Runge-Kutta 
method by the integration time step of 0.01 ms with double precision. Some results are 
examined by using the exponential method. We incorporate OU noises given by Eqs.(7)- 
(9) using the method of Fox, Gatland, Roy and Vemuri |57|. The initial conditions for 
the variables are given by 

V(t) = -65 mV, m{t) = 0.0526, h{t) = 0.600, n(t) = 0.313, for at t = 0, (10) 

which are the rest-state solution of a single HH neuron. We adopted parameters of 
V a = 30, V c = —50 mV, and r s = r n = 2 ms, with which the HH neuron is excitable 
receiving suprathreshold inputs. 



2.2 Wavelet Analysis 

2.2.1 Continuous Wavelet Transformation 

The continuous wavelet transformation (CWT) for a given regular function f(t) with a 
vanishing average is defined by 



c(a, b)=Jdt r a , b (t) f(t) =< i> a>b (t), f(t) >, (11) 

where a and b express the scale change and translation, respectively, and the star de- 
notes the complex conjugate. The wavelet function ip a ,b(t) is generated by dilation and 
translation of the mother wavelet ip{t) such as 

Mt) =| a I" 1 / 2 ^(^). (12) 

The parameters of a and b in Eqs.(ll) and (12) stand for the inverse of the frequency 
and the time, respectively. Then the CWT transforms the time-dependent function f(t) 
into the frequency- and time-dependent function c(a, b). The mother wavelet is a smooth 
function with good localization in both frequency and time spaces. A wavelet family 
given by Eq.(12) play a role of elementary function, representing the function f(t) as a 
superposition of wavelets ip a ,b{t)- 

The inverse of the wavelet transformation may be given by 

m=C?J ^Jdbc(a,b)Mt), (13) 

when the mother wavelet satisfies the following two conditions: 
(i) the admissibility condition given by 

< C < oo, (14) 

with 



C^= duj | | 2 / | u |, (15) 

J — oo 

where \l/(u;) is the Fourier transform of ip(t), and 



(ii) the zero mean of the mother wavelet: 

/°° 
df i[>(t) = *(0) = 0. (16) 
-oo 

Energy Distribution 

Parseval's theorem shows that the power spectrum is given by 

E tot = fdt | /(f) T= / ^ / db e(a, b) = J ' ^ E(a), (17) 
where the energy-density distribution e(a, b) and the a-dependent energy is given by 

e(a,b) = Cj 1 \c(a,b)\ 2 , (18) 

E(a) = J dbe(a,b). (19) 

Wavelet Entropy 

Equation (17) suggests that the probability of the power spectrum relevant to the 
states of a and b is given by 

p(a,b) = e(a,b)/E tat , (20) 

with the normalization condition: 

j*± Jdbp(a,b) = l. (21) 
We may define the wavelet entropy given by 

S= J dbs ( a > b )= j^ S ^ ( 22 ) 

where s(a,b) and Si (a) are given by 

s( a , b) = - p(a, b) log 2 [p(a, 6)], (23) 

Si(a) = J dbs{a,b), (24) 
Alternatively, we may define the wavelet entropy by 

S> = J da S 2 (a), (25) 

with 

S 2 (a) = - q(a) log 2 [g(a)], (26) 
where the probability q(a) is given by 

q(a) = a' 2 E(a)/E tot , (27) 
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satisfying the normalization condition: 

J da q(a) = 1, (28) 

the prefactor a~ 2 in Eq.(27) arising from the space factor relevant to the integration with 
respect to the variable a. Our numerical calculations to be presented in the following 
subsection shows that Si(a)/a 2 defined by Eq.(24) yields similar results as S 2 (a) defined 
by Eq.(26) 



2.2.2 Discrete Wavelet Transformation 

The discrete wavelet transformation (DWT) is defined for discrete values of a = 2 J and 
b = 2^k (j, k: integers) as 

c jk = c(2>, Vk) =< ^ fc (t), /(*) >, (29) 

with 

il> jk (t)=\2\- i ' 2 ij>(2- j t-k). (30) 
The ortho-normal condition for the wavelet functions is given by 

< ipjk(t), i/jj'k'(t) >= 8jj> 8 kk >, (31) 

which leads to the inverse DWT: 

/(0=E E Cik1> S k(t)- (32) 
j k 

In the multiresolution analysis (MRA) of the DWT, we introduce a scaling function 
(pit), which satisfies the recurrent relation with 2K masking coefficients, h k , given by 

2K-1 

<j>(t) = y/2 h k cf>(2t-k), (33) 

k=0 

with the normalization condition for <f>(t) given by 

J dt <j>{t) = I. (34) 
A family of wavelet functions are generated by 

2K-1 

iP(t) = E (-l) k h 2 K-i-k 0(2t - k). (35) 

k=0 

The scaling and wavelet functions satisfy the orthogonal relations: 

< 0(t), (f)(t - m) >= 8 m0 , < ijj(t),ip(t - m) >= 8 m0 , and < (f)(t),ip(t - m) >= 0. (36) 
The set of masking coefficients hj are chosen so as to satisfy the conditions shown above. 
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The simplest wavelet function for K = 1 is the Harr wavelet for which we get ho = 
h = l/y/2 and 

ip n (t) = 1 for < t < 1/2 (37) 
= -1 for 1/2 < t < 1, (38) 
= otherwise. (39) 

In the more sophisticated wavelets like the Daubechies wavelet, an additional condition 
given by 

[ dtt e ip{t) =0, for £ = 0,1, 2, 3 L-1 (40) 

is imposed for the smoothness of the wavelet function. Furthermore, in the Coiflet wavelet, 
for example, a similar smoothing condition is imposed also for the scaling function as 

J dt t l <f)(t) = 0, for £=1,2, 3 1! - 1 (41) 

In principle the expansion coefficients Cjf. in DWT may be calculated by using Eq.(29) 
for a given function f(t) and an adopted mother wavelet ip(t). This integration is, however, 
inconvenient, and in an actual fast wavelet transformation, the expansion coefficients are 
obtained by a matrix multiplication with the use of the iterative formulae given by the 
masking coefficients and expansion coefficients of the neighboring levels of indices, j and 
k 0. ' 

Energy Distribution 

Parseval's theorem for the DWT is given by 



E tot = jdt | f{t) | 2 = £ £ e jk = Y.E^ (42) 

3 k j 

where ejk, the energy distribution specified by j and k, and the j-dependent energy Ej 
are given by 

e?fc =| cjk | 2 , (43) 

Ej = Y,^- ( 44 ) 
k 

Wavelet Entropy 

We may define the wavelet entropy as 

S = Y,Sj, (45) 
j 

with 

Sj = ~J2 Pjk lo S2 Pjk, (46) 
k 

where the probability pjk is given by 

Pjk = e jk /E tot , (47) 
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satisfying the normalization condition: 



Vjk = 1- 

j k 

An alternative definition of the entropy is given by 0] |47 



with 



S'j = ~ qj log 2 qj, 



where qj expresses the probability given by 
with the normalization condition: 

E u = i- 



(48) 

(49) 

(50) 

(51) 
(52) 



A comparison between Sj and S'j will be numerically made in the following section. 

Hereafter, time, voltage, current, noise intensity (D), energy (E to t) and entropy (S) 
are expressed in units of ms, mV, /xA/cm 2 , //A 2 /cm 4 , mV 2 ■ ms and bits, respectively. 



3 Calculated Results 

Clustered spike trains to be analyzed are generated by using the HH neuron which receives 
input signals consisting of M (=1-5) pulses with the ISI of T im = 25 ms (m=l-4). Pulse 
clusters are separated by the longer interval of 100 ms. We choose this ISI value of 25 
ms because 7-wave spikes with the average frequency of about 40 Hz are ubiquitous in 
brain [p7|-[p9l. The time dependence of U- u input potential I (= I ps + J n ) and membrane 



potential V is shown in Fig. 1. Spikes fire after an injection of input pulses with a delay 
of about 2 ~ 3 ms. 

Generated spike trains have been analyzed with the use of the WT. One of the advan- 
tages of the WT over the FT is that we can choose a proper mother wavelet depending 
on the shape of signals to be examined. Among many candidates of mother wavelets, we 
have adopted the Coiflet because its shape is similar to that of spikes. Compromising the 
accuracy and the computation effort, we have decided to adopt the Coiflet of order 3 as 
the mother wavelet for our analysis. 

We analyze the time- dependent membrane potential V(t) by setting the signal to be 
f{t) = V(t)- < V(t) > in Eqs.(ll) and (29), < V(t) > being the time average of V(t). 
Spike signals V(t) are assumed to be given at the sampled t values with a sampling time 
of At = 1 ms, and are analyzed by both the CWT and DWT with the use of MATLAB 
wavelet tool box. 

The analyzed spike train is depicted in Fig. 2(a) and the pattern of the expansion 
coefficients c(a,b) obtained by the CWT is shown in the (a, b) space of Fig. 2(b), where 
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the medium dark tone denotes the zero level and the black (white) expresses the nega- 
tive (positive) value. Note that 1/a and b physically stand for the frequency and time, 
respectively. The pattern of the energy- density distribution e(a, b) [=| c(a, b) | 2 by setting 

= 1 in Eq.(18) hereafter] is shown in Fig. 2(c), where black and white denote the 
minimum (zero) and maximum levels, respectively. Patterns near the both sides along 
vertical axis in Fig. 2(b) and (c) are due to the edge effect and should not be taken seri- 
ously. Figure 2(c) shows that for a single (M = 1) spike signal, the energy- density pattern 
shows a cone which is progressively broader at larger a. For spike signals with multiple M, 
cones of spikes are separable for small a whereas they merge for larger a, yielding intrigue 
structure in the pattern. Figure 3 shows the b dependence of the expansion coefficients 
c(a, b) for a = 2, 4, 8, 16, 32 and 64. We notice that the CWT coefficients for smaller a 
oscillates with shorter period. The magnitudes of c(a, b) are largest for a ~ 8 — 16 n Fig. 
3: note vertical scales for a— 2 and 64 to be different from those for the other a values. 

In what follows, we pay our attention to the M = 3 spike in order to perform more 
detailed analysis. Figure 4(a) shows the M = 3 spike signal to be analyzed. The energy- 
density distribution analyzed by the CWT is depicted in Fig. 4(b) . At small a each pulse 
yields a characteristic cone pattern. At 10 ~ a ~ 20, cones arising the consecutive pulses 
separated by T Q = 25 msec, begin to merge. At larger a ~ 40 cones originating from the 
first and the last pulses, which are separated by 2T D = 50 msec, begins to correlate. The 
pattern has structures near both end sides of vertical axis because of the edge effect. 

Figure 5 expresses the b dependence of the expansion coefficients of c(a, b) for typical 
a values calculated in the CWT. As noticed in Fig. 3, the magnitudes of c(a, b) is most 
significant for a ~ 16. This is more clearly seen in the a-dependence of the energy- density 
distribution E(a) shown in Fig. 6(a) which has a peak at a ~ 20. However, when taking 
into account the space factor a" 2 in the integration with respect to the variable a [see 
Eq.(17)], we note that E(a)/a 2 has a peak at much a lower value of a = 4. 

The a-dependent average \p(a)] and root-mean-square value [a (a)} of CWT coefficients 
are plotted in Fig. 7(a) and (b), respectively. A rise in /i(a) at a > 30 is due to the edge 
effect. We note that a (a) has a peak at a ~ 20. When the space factor a~ 2 is included, 
the peak position of a(a)/a 2 moves to the lower value of a = 4. 

So far we have employed the CWT. When we apply the DWT to the signal shown in 
Fig. 4(a), we get the j— and A;-dependent energy distribution, e jft , which is plotted in 
Fig. 4(c). Because the scale (a = 2 j ) and translation (b = 2 j k) parameters are discrete 
in the DWT, the energy distribution, e^, defined by Eq.(43) are depicted by the blocks 
in the (a, b) space. The DWT decomposition of the spike signal: V = Z)j=i a j + ^6, is 
shown in Fig. 8. A contribution from j = 2 (a = 4) to V seems to be predominant. This 
is supported by the calculated j-dependent energy Ej plotted in Fig. 6(c), which has the 
maximum at j = 2 (a = 4). The j dependence of the average (fj,j) and RMS value (aj) 
of the DWT coefficients are plotted in Fig. 7(c) and (d), respectively. The maximum in 
aj is realized at 2? = 8 (j = 3), which should be compared with the maximum at a = 4 
in a(a)/a 2 . These comparisons show that the a dependence of E(a)/a 2 [a(a)/a 2 ] in the 
CWT are in fairly good agreement with the j dependence of Ej [aj] in the DWT. 

Next we discuss the wavelet entropy calculated in terms of the WT expansion coeffi- 
cients. Figure 9(a) shows the contour map depicting the entropy density s(a, b) defined by 
Eq.(23) with the CWT. Figure 9(b) shows the a dependence of s(a,b) for typical values 



11 



of a = 103, 109 and 115, which are indicated by vertical, dashed lines in Fig. 9(a). s(a, b) 
has peaks at a ~ 10, a ~ 15 and a ~ 20 for b = 103, b = 109 and b = 115, respectively. 
Figure 9(c) expresses the b dependence of s(a, b) for b = 10, 20 and 30 ms, which are shown 
by horizontal, dashed lines in Fig. 9(a). Although s(a,b) for a = 10 has the maximum 
peaks at the firing times t on (=103, 128 and 153 ms), it has also side-peaks before and 
after t on . For a = 20 and 30, s(a, b) has peaks not only at t on but also between the firing 
times. The a dependence of the CWT entropy Si (a) defined by Eq. (24) is plotted in Fig. 
6(b), which has a peak at a = 20, just as E(a) shown in Fig. 6(a). We note, however, 
that Si(a)/a 2 has a peak at a lower value of a = 3, which may be compared with the peak 
position at a = 4 of the entropy S 2 (a) defined by Eq.(26). 

The DWT entropies, Sj and Sj, defined by Eqs.(46) and (50), are shown in Fig. 6(d). 
Although Sj is about two times larger than Sj, both entropies have a similar j dependence 
with the maximum at a = 4 (j = 2). 

3.1 Effects of Variation in ISI 

We will discuss in this subsection, how the variation in ISI of the spike signal modifies 
the result of the WT. Firstly we equally change both T ol and T o2 as T ol = T o2 = T a . 
Figure 10(a), (b), (c) and (d) show the analyzed spike signals (upper frames) and the 
energy-density patterns (lower frames) for T Q = 15, 25, 35 and 50 ms, respectively. We 
note that the size of patterns becomes large almost in proportion to an increase in T Q . 
Figure 11(a) and (b) show the a dependence of the energy-density distribution E(a) and 
the entropy S(a) in the CWT, respectively. The peak position moves to larger a for the 
spike signal with larger T Q , as expected. The relevant j dependence of Ej and Sj in the 
DWT is shown in Fig. 11(c) and (d). When T Q is smaller, Ej and Sj have large magnitude 
for smaller j, which is consistent with the result obtained in the CWT. 
The FT for a signal fit) is defined by 

F{oo)= fdte-™f(t). (53) 

When f(t) is given by a delta- function- type function: 

/(*) = E *(f-fon)> (54) 
n=l 

which is a simplification of our M = 3 HH spike with firing times of t on , its FT spectrum 
is given by 

|FH| 2 =2 ]T cos(27r/// on ) + 3. (55) 

n=l 

In Eq.(55) / = u/2rr is the frequency and the fundamental frequencies f on (n — 1 — 3) 
are given by f ol = l/T ol , f o2 = l/T o2 and / o3 = 1/ (T ol + T o2 ). 

Figures 12(a), (b) (c) and (d) show the FT frequency spectra of the HH spike signals 
for T = 15, 25, 35 and 50 ms, respectively. When T Q = 25, for example, the fundamental 
frequencies are / ol = f o2 = 40 and / O 3=20 Hz, the former being coincidentally the second 
harmonics of the latter. The dashed curve in Fig. 12(b) shows the FT spectra given by 
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Eq.(55) for a delta-function- type signal given by Eq.(54) with T G = 25 ms, which is in 
good agreement with the FT spectra of the HH spike at / < 150 Hz. 

Next we investigate the effect of ISI on the WT result by changing T a \ but keeping 
T ol + T o2 = 50 ms. The analyzed spike signals (upper frames) and the pattern of the 
energy-density distribution (lower frames) for T ol = 10, 15, 20 and 25 ms are shown in 
Fig. 13(a), (b), (c) and (d), respectively. Note that the average ISIs are the same (25 
ms) but their RMS are different: (< (T Q 2 n - < T on > 2 ) >) 1 ' 2 = 15, 10, 5 and ms for T ol 
= 10, 15, 20 and 25 ms, respectively. For the case of T ol (= T o2 ) = 25 ms, the calculated 
pattern in Fig. 13(d) is nearly symmetric with respect to the axis of b = 128 ms, at which 
the second spike fires. With introducing a small asymmetry in ISI: T ol = 20 ( T o2 = 30) 
ms, the pattern becomes asymmetric, as shown in Fig. 13(c). This is more clearly seen 
in Fig. 13(a) for the case of T Q \ = 10 ( T o2 = 40) ms: smaller T i and larger T o2 yield 
smaller and larger patterns, respectively. Figure 14(a) and (b) show the a dependence of 
E(a) and S(a), respectively. In the case of T ol = 25 ms, E(a) has a single maximum at 
a ~ 20. On the contrary, in the case of T ol = 10 ms, E(a) and S(a) have double maxima: 
the peak at smaller a arises from the contribution due to a smaller T ol and the peak at 
larger a from a larger T o2 . These features are consistent with the results obtained in the 
DWT, which are shown in Figs. 14(c) and 14(d). Figures 13 and 14 clearly show that the 
response of the HH neuron cannot described only by the rate of the input pulses [5f| , in 



agreement with experiments |)9| . 

Figure 15 shows the FT spectra for spike signals for T ol = 10, 15, 20 and 25 ms. In the 
case of T i = 15 ms, for example, the fundamental frequencies are 66.7, 28.6 and 20 Hz. 
Although a glimpse of Fig. 15(b) suggests that fundamental frequencies are 60 and 80 
Hz, it is not true. The dashed curve in Fig. 15(b) shows the FT spectra given by Eq.(55) 
for a delta-function signal given by Eq.(54) with / ol = 66.7, f o2 = 28.6 and f o3 = 20 Hz, 
which is in good agreement with the FT spectra of HH spikes at / < 150 Hz. 



3.2 Effects of Noises 

We will study, in this subsection, the effect of OU noises given by Eqs.(7)-(9). Figure 
16(a), (b), (c) and (d) show the spike signals (upper frames) and the CWT patterns of 
the energy- density distribution (lower frames) for D — 0, 1, 2 and 3, respectively. The 
result without noises (D = 0) has been shown previously [e.g. Fig. 4(b)]. When noises 
are introduced, spike signals lead to fine structures in the CWT patterns. For a large 
D = 3, applied noises trigger a spurious spike at t — 203 ms. 

Figure 17(a) and (b) show E(a) and Si (a) for D = 0, 1, 2 and 3 calculated by the 
CWT coefficients. Applied noises provide extra energy contributions in a fairly wide range 
of a value. Particularly for D = 3, E(a) and S(a) are much increased by a spuriously 
triggered spike at t — 203 ms. 

Similar results are obtained in the j-dependent energy (Ej) and entropy (Sj) calculated 
by the DWT, which are plotted in Fig. 17(c) and (d). 

Spike signals with noises are analyzed also by the FT, whose frequency spectra for 
D = 0, 1, 2 and 3 are plotted in Fig. 18(a), (b), (c) and (d), respectively. It is interesting 
that peaks of the fundamental frequency of 40 Hz and its harmonics are sharpened by 
noises. 
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4 Conclusion and Discussion 

One of the advantages of the DWT is that we can easily separate noise components from 
a given signal. As the first example, we discuss the denoising by using the DWT.f60|j- [|63|j 
The key point in the denoising is how to choose which wavelet coefficients are correlated 
with the signal and which ones with noises. The simple denoising is to neglect some DWT 
expansion coefficients when reproducing the signal by the inverse wavelet transformation. 
In getting the inverse WT by 

/iW = EE4W). ( 56 ) 

3 k 

we may assume that the components for a < a c in the (a, b) plane arise from noises to set 
the denoising coefficients as 

C S = °jki for 3 > 3c (a > a c ) 

= 0, otherwise (Method I) (57) 

where j c = log 2 a c is the critical j value. A demonstration of denoising of spike signals 
with M = 3 and D = 3 is given in Fig. 19. Figure 19(a) expresses the original output 
spike, and Fig. 19(b) shows the inverse signals by using the method I [eq.(4.2)] with 
j c = 2. For a comparison, we show the relevant results for the postsynaptic input I{. the 
original signal in Fig. 19(d) and the inverse signal with denoising for j c = 2 in Fig. 19(e). 

Since the DWT transforms the original signal to the two-dimensional (a, b) plane, 
it has much freedom than the FT, which makes it possible for us to adopt the more 
sophisticated denoising. We may assume that the components for b < b^ or b > by at 
a < a c in the (a,b) plane are noises to set the denoising WT coefficients as 

c fk = c jk, for 3 < 3c or k L < k < k v 

= 0, otherwise (Method II) (58) 

where j c = log 2 a c , and k^ = 2^ and k\j = b\j 2~i are lower and upper critical k values. 
Figures 19(c) and 19 (f) show the donoising signals of the output and postsynaptic signals, 
respectively, by using the method II [Eq.(4.3)] with 6 L = 40, &u = 220 and j c = 4. The 
denoising method II given by Eq. (4.3) is better than the method I given by eq.(4.2) 
because the magnitude of a spurious spike at t ~ 200 ms is much reduced while the signal 
at 100 < t < 160 ms is better reproduced in the method II than in the method I. 

As the second example using the DWT, we discuss the signal-to-noise ratio (SNR). 
From the above consideration, we may define the signal component A s and the noise 
component A Q by 

^ = EEI4I 2 > ( 59 ) 

3 k 

4 = EDI^ \ 2 ~\c%\ 2 ) (60) 

3 k 

which lead the SNR r\: 

r l = \og w (A s /A n ). (61) 
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The solid curves in Fig. 20(a) and Fig. 20(b) shows the -D-dependent SNR of output 
(V) and input signals (ij), respectively, with M = 3 calculated by the denoising method 
I given by eq.(4.2) with j c = 2. It might be strange that SNR is finite even for no noises 
(D = 0). This is because signals for D = already included in the j = 1 component 
are regarded as noises in the denoising method I though its magnitude is small. On the 
contrary, the dashed curves in Fig. 20(a) and 20(b) show the SNR as a function of D 
calculated by the denoisng method II given by eq.(4.3) with b L = 40, bjj = 220 and 
j c = 4. Because of the condition imposed for the a variable in the method II, the noise 
contribution is much reduced for small D, yielding a large SNR. With increasing the value 
of D, SNR of the input and output signals decreases as expected. 

Finally, we apply the CWT to a signal which is a simplified one of the original HH 
neuron spike, in order to get some insight how the detailed shape of spikes is relevant to 
information transmission. As its candidate, we adopt a square pulse with the width of 1 
ms, the delta- function- type signal. Figures 21(a) and (b) show the analyzed signal and its 
energy-density distribution in the CWT, respectively. Comparing Fig. 21(b) with Fig. 4 
(b) for the original HH spike, we note that except for very small value of a(~ 4), the CWT 
pattern for the delta-function-type signal is almost the same as that of the original HH 
spike. This suggests that information is predominantly carried by firing times of spikes, 
but not by their detailed structure. This agrees with the conventional wisdom in the 
neuroscience community. 

To summarize, we have analyzed transient spike signals with the use of both the CWT 
and DWT. The WT has been shown to be more useful than the FT for an analysis of 
transient signals like spikes. The WT now has a wide-range of applications including the 
information compression like MPEG and JPEG. It might be possible that real neural 
systems adopt the WT-like technique for their efficient information transaction. 
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Figure Captions 

Fig. 1 The time course of (a) the input potential U h (b) the postsynaptic potential / 
(= -fps + ^n) and (c) the membrane potential V . 

Fig. 2 (a) The analyzed signal, (b) the pattern of the CWT coefficients c(a,b), and (c) 
that of the energy-density distribution e(a, b) =| c(a, b) \ 2 . Scales of colors from minimum 
to maximum in (b) and (c) are shown. 

Fig. 3 The b dependence of the CWT coefficients c(a, b) for specified a values. 

Fig. 4 (a) The M = 3 spike signal, (b) the profile of energy distribution e(a, b) analyzed 
by the CWT, and (c) that analyzed by the DWT. 

Fig. 5 The b dependence of the CWT coefficients c(a, b) for specified a values. 

Fig. 6 (a) The a-dependent energy [E(a) and E(a)/a 2 ] and (b) entropy [S(a) and S(a)/a 2 ] 
in CWT. (c) The j-dependent energy (Ej) and (d) entropy (Sj, £J) in DWT, solid lines 
being drawn only for a guide of eye (see text). 

Fig. 7 (a) The a-dependent average of c(a,b) [//(a) and /i(a)/a 2 ], and (b) their RMS 
[cr(a) and a /a 2 ] in CWT. (c) The j-dependent average of Cjk (fJ-j), and (d) their RMS (cr,-) 
in DWT, solid lines being drawn only for a guide of eye (see text). 

Fig. 8 The DWT decomposition of the signal to various components; V = X^=i dj + o,q. 

Fig. 9 (a) The contour map of the entropy distribution s(a, b). (b) The a dependence of 
s(a, b) for b = 103 (solid curve), 109 (dot-dashed curve) and 115 ms (dashed curve), whose 
positions are indicated by vertical, dashed lines in (a), (c) The b dependence of s(a, b) for 
a = 10 (solid curve), 20 (dot-dashed curve) and 30 (dashed curve), whose positions are 
indicated by horizontal, dashed lines in (a). 

Fig. 10 The time dependence of spikes (upper frames) and the energy-density patterns 
(lower frames) in the CWT for (a) T G = 15, (b) 25, (c) 35 and (d) 50 ms. 

Fig. 11 (a) The a dependence of E(a) and (b) the entropy Si (a) in the CWT for T Q = 
15, 25, 35 and 50 ms. (c) The j dependence of Ej and (d) the entropy Sj in the DWT for 
T = 15, 25, 35 and 50 ms. 

Fig. 12 The FT spectra for (a) T D = 15, (b) 25, (c) 35 and (d) 50 ms. The dashed curve 
in (b) expresses the FT spectra given by Eq.(55) for a delta-function-type signal given by 
Eq.(54) with T G = 25 ms. 

Fig. 13 The time dependence of spikes (upper frames) and the energy-density patterns 
(lower frames) in the CWT for (a) T oi = 10, (b) 15, (c) 20 and (d) 25 ms. 
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Fig. 14 (a) The a dependence of E(a) and (b) the entropy Si (a) in the CWT for T \ = 
10, 15, 20 and 25 ms. (c) The j dependence of Ej and (d) the entropy Sj in the DWT 
for T i = 10, 15, 20 and 25 ms. 

Fig. 15 The FT spectra for (a) T ol = 10, (b) 15, (c) 20 and (d) 25 ms. The dashed curve 
in (b) expresses the FT spectra given by Eq.(55) for a delta-function-type signal given by 
Eq.(54) with T ol = 15 ms. 

Fig. 16 The time dependence of spikes (upper frames) and the energy-density patterns 
(lower frames) in the CWT for (a) D = 0, (b) 1, (c) 2, and (d) 3. 

Fig. 17 (a) The a dependence of E(a) and (b) the entropy Si (a) in the CWT for D = 
0, 1, 2 and 3. (c) The j dependence of Ej and (d) the entropy Sj in the DWT for D = 
0, 1, 2 and 3. 

Fig. 18 The FT spectra for (a) D = 0, (b) 1, (c) 2, and (d) 3. 

Fig. 19 (a) The original output spike V with D = 3, (b) its inverse signal by the denoising 
method I [eq.(4.2)] with j c = 2, and (c) that by the denoising method II [eq.(4.3)] with 
&l = 40, b\j = 220 and j c = 4. (d) The original input signal with D — 3 (e) its inverse 
signal by the denoising method I [eq.(4.2)] with j c = 2, and (f) that by the denoising 
method II [eq.(4.3)] with b L = 40, b v = 220 and j c = 4. 

Fig. 20 The D dependence of SNR of (a) the output spike signal V and (b) the input 
presynaptic signal the solid curves show the results by using the denoising method I 
[eq.(4.2)] with j c = 2, and the dashed curves express the results by using the denoising 
method II [eq.(4.3)] with b L = 40, b v = 220 and j c = 4 (see text). 

Fig. 21 (a) The analyzed delta-function-type signal and (b) the energy-density pattern 
in the CWT (see text). 
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